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[i]  Difficulties  associated  with  identifying  the  dense  nonaqueous  phase  liquid  (DNAPL) 
source  zone  architecture  at  the  field  scale,  combined  with  the  computational  costs  of 
field-scale  DNAPL  dissolution  simulations,  have  motivated  the  development  of  a  number 
of  simplified  models  that  rely  upon  upscaled  (i.e.,  domain-averaged)  mass  transfer 
coefficients  to  approximate  field-scale  dissolution  processes.  While  conceptually 
attractive,  these  upscaled  models  have  yet  to  be  fully  evaluated  for  prediction  of  mass 
recovery  from  a  range  of  nonuniform,  three-dimensional  DNAPL  source  zones.  This  study 
compares  upscaled  model  predictions  of  flux-weighted  downstream  concentrations  and 
source  longevity  to  predictions  derived  from  three-dimensional  multiphase  numerical 
simulation  of  tetrachloroethene  (PCE)-NAPL  dissolution  for  realizations  of  a  statistically 
homogeneous,  nonuniform  aquifer.  Although  the  functional  forms  of  the  upscaled  models 
are  generally  shown  to  be  mathematically  equivalent,  upscaled  model  flux-weighted 
concentration  predictions  varied  by  over  one  order  of  magnitude,  with  variations  attributed 
to  the  dependence  of  the  upscaled  model  parameters  on  the  specific  source  zone  scenario 
used  for  model  calibration.  Replacement  of  upscaled  model  calibration  parameters  with 
source  zone  parameters  that  can  be  obtained  from  site  characterization  information 
(specifically,  the  initial  flux-weighted  concentration  and  source  zone  ganglia-to-pool 
(GTP)  mass  ratio)  reduced  the  root-mean-square  error  between  upscaled  and  numerical 
model  predictions  by  approximately  80%.  Application  of  this  modified  model  to  a  range 
of  source  zone  scenarios  (0.4  <  GTP  <  oo)  demonstrates  the  efficacy  of  the  model  for 
use  as  a  screening  tool  to  relate  DNAPL  mass  removal  and  flux-weighted  concentrations 
when  mass  removal  is  less  than  80%. 
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1.  Introduction 

[2]  The  presence  of  dense  nonaqueous  phase  liquid 
(DNAPL)  at  a  site  is  a  critical  factor  in  the  design  and 
efficacy  of  any  site  remediation  strategy  [Mackay  and 
Cherry,  1989].  When  released  into  an  aquifer,  a  DNAPL 
can  become  entrapped  as  discontinuous  ganglia  [C/.5. 
Environmental  Protection  Agency  ( USEPA ),  1990;  Mercer 
and  Cohen,  1990]  and  accumulate  in  higher  saturation  pools 
at  textural  interfaces  within  the  medium  [Kueper  et  al, 
1993],  serving  as  a  long-term  source  of  dissolved  phase 
contamination.  The  region  of  a  contaminated  aquifer  con¬ 
taining  DNAPL  pools  and  ganglia  is  frequently  referred  to 
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as  a  source  zone  [. National  Research  Council  ( NRC ),  2005]. 
Recently,  research  has  focused  on  quantifying  the  potential 
benefits  of  DNAPL  source  zone  mass  removal,  including 
reductions  in  down-gradient  contaminant  concentrations, 
mass  flux,  and  DNAPL  source  longevity  [ Berglund ,  1997; 
Dekker  and  Abriola,  2000b;  Sale  and  McWItorter,  2001; 
Rao  et  al.,  2001;  Rao  and  Jawitz,  2003;  USEPA, 
2003;  Lemke  et  al.,  2004b;  Parker  and  Park,  2004;  Soga 
et  al.,  2004;  Wood  et  al.,  2005;  Jawitz  et  al.,  2005;  Fure  et 
al.,  2006],  A  number  of  researchers  have  developed  screen¬ 
ing  models  that  relate  reductions  in  down-gradient  contam¬ 
inant  concentration  and  source  longevity  to  the  level  of 
mass  removal  [Dekker,  1996;  Parker  and  Park,  2004;  Zhu 
and  Sykes,  2004;  Falta  et  al.,  2005a;  Park  and  Parker, 
2005;  Jawitz  et  al.,  2005;  Fure  et  al.,  2006],  These  models 
are  typically  simplified  parametric  or  analytical  expressions 
of  reduced  dimensionality.  The  development  of  these  mod¬ 
els  has  been  motivated  by  difficulties  associated  with 
applying  more  comprehensive  modeling  tools,  including 
description  of  DNAPL  source  zone  architecture  at  the  field 
scale  and  computational  costs  of  field-scale  dissolution 
simulations.  The  general  utility  of  these  screening  models, 
however,  is  unclear  given  that  they  cannot  be  used  in  a 
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predictive  sense  due  to  the  need  to  calibrate  model  param¬ 
eters  to  specific  site  conditions. 

[3]  Simplified  screening  models  may  be  divided  into  two 
categories:  (1)  stochastic -advective  models,  which  relate 
source  zone  remediation  performance  to  joint  spatial  vari¬ 
ability  of  both  groundwater  flow  and  NAPL  content 
[Berglund,  1997;  Jawitz  et  al,  2003,  2005;  Enfield  et  al., 
2005;  Wood  et  al,  2005;  Fure  et  al.,  2006],  and  (2)  upscaled 
models  (sometimes  referred  to  as  macroscopic,  effective,  or 
field  scale),  which  incorporate  the  effects  of  spatial  varia¬ 
tions  in  DNAPL  saturations  and  flow  bypassing  on  source 
zone  mass  discharge  through  the  use  of  an  upscaled  (i.e., 
domain-averaged)  mass  transfer  coefficient  [ Parker  and 
Park,  2004;  Zhu  and  Sykes,  2004;  Falta  et  al.,  2005a;  Park 
and  Parker,  2005;  Basu  et  al.,  2006].  Stochastic-advective 
models  have  been  used  to  approximate  the  response  in 
contaminant  mass  flux  due  to  DNAPL  source  depletion 
using  trajectory-averaged  model  parameters  [ Wood  et  al, 
2005;  Fure  et  al.,  2006].  These  trajectory-averaged  model 
parameters  are  flow  field  specific  and  are  typically  obtained 
by  interrogating  the  swept  volume  with  a  tracer  or  suite  of 
tracers  [Jawitz  et  al.,  2003;  Enfield  et  al.,  2005].  In  contrast, 
upscaled  modeling  methods  have  been  employed  to  quan¬ 
tify  the  relationship  between  DNAPL  mass  and  source 
strength  [ Falta  et  al.,  2005a,  2005b;  Basu  et  al.,  2006]. 
Application  of  these  upscaled  models  does  not  require 
imposition  of  the  remedial  flow  field  and  may  thus  prove 
better  suited  for  incorporating  the  information  gained  during 
a  typical  site  characterization  effort  into  a  remedial  design. 
The  present  work  is  restricted  to  consideration  of  upscaled 
modeling  methods  [Parker  and  Park,  2004;  Zhu  and  Sykes, 
2004;  Falta  et  al.,  2005a;  Park  and  Parker,  2005]. 

[4]  Upscaled  models  typically  assume  that  the  flow  field 
is  steady  and  dissolution  mass  transfer  can  be  approximated 
by  a  linearized  Fick’s  law  expression  incorporating  an 
upscaled  (i.e.,  domain-averaged)  mass  transfer  coefficient. 
Upscaled  mass  transfer  coefficients  have  been  estimated 
using  explicit  descriptions  of  pore  networks  and  NAPL 
saturation  distribution  [Jia  et  al.,  1999;  Field  and  Celia, 
2001;  Sale  and  McWhorter,  2001],  or  by  fitting  solutions  of 
a  simplified  one-dimensional  component  mass  balance 
equation  to  experimental  [Saba  and  Illangasekare,  2000; 
Schaerlaekens  and  Feyen,  2004]  or  numerically  generated 
results  [Dekker,  1996;  Parker  and  Park,  2004;  Zhu  and 
Sykes,  2004;  Park  and  Parker,  2005],  In  the  latter  cases,  the 
upscaled  mass  transfer  coefficient  has  generally  been  cor¬ 
related  with  several  simulation  parameters  including 
groundwater  velocity,  mean  hydraulic  conductivity,  and 
the  extent  of  mass  removal  using  best  fit  correlation 
parameters.  These  upscaled  mass  transfer  correlations  are 
similar  in  form  to  correlations  that  have  been  derived  in 
other  studies  to  describe  local-scale  dissolution  under  a 
variety  of  conditions,  including  dissolution  of  residual 
NAPL  ganglia  [e.g.,  Miller  et  al.,  1990;  Geller  and  Hunt, 
1993;  Powers  et  al.,  1994;  Imhoff  et  al,  1994],  dissolution 
from  high  initial  saturation  ganglia  regions  [e.g.,  Nambi  and 
Powers,  2000,  2003],  interphase  mass  exchange  from 
DNAPL  pools  [e.g.,  Kim  and  Chrysikopoulos,  1999; 
Chrysikopoulos  and  Kim,  2000],  and  dissolution  from  a 
source  emplaced  in  the  field  [Frind  et  al.,  1999;  Broholm  et 
al.,  1999,  2005;  Rivett  and  Feenstra,  2004],  Unlike  the 
local-scale  mass  transfer  correlations,  which  have  demon¬ 


strated  applicability  for  a  range  of  porous  media,  flow  and 
entrapment  conditions  [e.g.,  Powers  et  al.,  1994;  Imhoff  et 
al.,  1994],  existing  upscaled  mass  transfer  correlations  tend 
to  be  site  specific,  valid  only  for  sites  with  conditions  (i.e., 
source  zone  architectures)  consistent  with  those  used  to 
parameterize  the  upscaled  mass  transfer  correlation  [Parker 
and  Park,  2004]. 

[5]  Despite  their  simplicity  and  conceptual  attractiveness, 
the  utility  of  the  upscaled  model  as  an  alternative  to 
computationally  intensive  numerical  simulations  is  limited 
unless  upscaled  mass  transfer  coefficient  correlations  can  be 
developed  that  apply  to  a  broad  range  of  sites.  The  objec¬ 
tives  of  this  work  are  (1)  to  demonstrate  the  predictive 
limitations  of  existing  upscaled  mass  transfer  models, 
through  comparison  with  three-dimensional  numerical  sim¬ 
ulations  of  dissolution  from  representative  DNAPL  source 
zones,  (2)  to  develop  an  improved  formulation  that  facili¬ 
tates  independent  estimation  of  model  parameters  over  a 
wide  range  of  source  zone  scenarios,  and  (3)  to  evaluate  the 
predictive  capability  of  this  improved  upscaled  mass  trans¬ 
fer  model.  The  source  zone  scenario  selected  for  these 
evaluations  was  based  on  the  Bachman  Road  site,  which 
was  extensively  characterized  as  part  of  a  pilot-scale  sur¬ 
factant  enhanced  aquifer  remediation  (SEAR)  demonstra¬ 
tion  [Abriola  et  al,  2005;  Ramsburg  et  al.,  2005],  This  site 
is  representative  of  many  small-scale  DNAPL  source  zones 
located  throughout  the  United  States  [Abriola  et  al.,  2005]. 

2.  Methodology 

[6]  The  objectives  of  this  work  require  the  use  of  both 
numerical  and  analytical  methods  to  model  mass  discharge 
from  DNAPL  source  zones.  Comprehensive  reviews  of 
compositional  multiphase  numerical  models  can  be  found 
in  the  literature  [Abriola,  1989;  Miller  et  al,  1998;  Adeel  et 
al.,  2001],  Likewise,  analytical  models  describing  contam¬ 
inant  fate  and  transport  are  discussed  in  most  groundwater 
textbooks  [e.g.,  Bear,  1979;  Fetter,  1993],  Therefore  this 
section  includes  only  a  brief  review  of  the  equations 
governing  multiphase,  multicomponent  flow.  Simplifying 
assumptions,  equation  solution  techniques,  and  site  condi¬ 
tions  employed  in  the  simulations  are  presented  for  both  the 
numerical  and  analytical  models. 

2.1.  Governing  Equations 

[7]  Simulation  of  a  DNAPL  spill  and  subsequent  con¬ 
taminant  mass  discharge  in  the  saturated  zone  requires  the 
solution  of  both  phase  and  component  mass  balance  equa¬ 
tions.  Phase  mass  balance  equations  describe  bulk  fluid 
phase  flow  using  a  modified  form  of  Darcy’s  equation  and 
constitutive  relationships  for  saturation,  relative  permeabil¬ 
ity,  and  capillary  pressure  [see  Abriola,  1989],  Within  each 
bulk  fluid  phase  (aqueous,  a  =  a;  organic,  a  =  n),  the 
spatial-temporal  distribution  of  component  i  is  computed 
using  a  component  mass  balance  equation  written  in  terms 
of  the  mass  concentration  of  component  i  in  the  a  phase 

(c?y. 

*4  C*«c?)  +  <^V  ■  sa(C?Va  -  D‘a  ■  VC?)  =  (!) 

u  p 

<j>  is  the  matrix  porosity,  sa  is  the  a  phase  saturation,  D'a  is 
the  three-dimensional  hydrodynamic  dispersion  tensor  for 
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component  i  in  phase  a  [Bear,  1972],  V'  is  the  a  phase  pore 
velocity  computed  using  a  modified  form  of  Darcy’s  law 
[Abriola,  1989],  and  Eapi  is  the  interphase  mass  exchange 
of  component  i  from  the  3  to  the  a  phase  [Weber  and 
DiGiano,  1996]: 

EaPi  =  K*[Si{C?eq-C“)  (2) 

where  c?e<!  is  the  equilibrium  solubility  of  component  i  in  the 
a  phase,  and  nal 3,:  is  a  lumped  mass  transfer  coefficient, 
which  quantifies  the  rate  of  mass  transfer  of  component  i 
from  the  (3  to  the  a  phase.  This  lumped  coefficient  («aig;)  is 
generally  represented  as  a  correlation  between  a  Sherwood 
number  (NSh  —  napid\olDm),  normalized  DNAPL  saturation 
Reynolds  number  ( Nrb  =  Va pad50/ pf),  and  Schmidt 
number  (Nsc  =  p JpaDm ),  where  d50  is  the  mean  grain 
diameter,  and  Dm  is  the  aqueous  phase  molecular  diffusion 
coefficient  [e.g.,  Miller  et  al.,  1990;  Imhoff  et  al.,  1994; 
Powers  et  al,  1994],  Incorporation  of  the  normalized 
saturation  into  the  lumped  mass  transfer  coefficient 
accounts  for  the  transient  nature  of  the  lumped  mass 
transfer  coefficient  resulting  from  reductions  in  interfacial 
area  as  the  DNAPL  ganglia  dissolve. 

[8]  Solution  of  the  phase  and  component  mass  balance 
equations  permits  the  quantification  of  several  commonly 
employed  source  zone  management  metrics  which  may  be 
used  to  assess  the  performance  of  the  upscaled  model: 
( 1 )  flux-weighted  concentration  at  a  down-gradient  boundary 
or  point  of  compliance,  (2)  contaminant  mass  flux  across  a 
down-gradient  boundary,  and  (3)  source  longevity.  These 
metrics  have  been  used  extensively  to  quantify  the  benefits  of 
partial  mass  removal  from  DNAPL  source  zones  [Dekker  and 
Abriola,  2000b;  Einarson  and  Mackay,  2001;  Sale  and 
McWhorter,  2001;  McWhorter  and  Sale,  2003;  Rao  et  al., 
2001;  Saenton  et  al.,  2002;  Rao  and  Jawitz,  2003;  Lemke  et 
al.,  2004b;  Soga  et  al.,  2004;  Broholm  et  al.,  2005;  Jawitz  et 
al.,  2005;  Fure  et  al.,  2006;  Basu  et  al.,  2006],  and  have 
recently  been  proposed  as  an  alternative  to  more  traditional 
remediation  metrics,  such  as  maximum  contaminant  concen¬ 
tration  [Stroo  et  al.,  2003;  USEPA,  2003;  NRC,  2005]. 

2.2.  Numerical  Simulations 

[9]  Source  zone  distributions  have  generally  been  char¬ 
acterized  using  domain-averaged  NAPL  saturations,  sto- 
chastic-advective  trajectory  averaged  NAPL  saturations 
[e.g.,  Jawitz  et  al.,  2005],  or  ganglia-to-pool  (GTP)  mass 
ratios  [Lemke  et  al.,  2004b;  Christ  et  al.,  2003;  Lemke  and 
Abriola,  2006].  While  each  of  these  metrics  has  utility,  the 
GTP  mass  ratio  is  selected  for  analysis  herein  due  to  its 
potential  to  link  mass  discharge  strength  to  source  zone 
architecture  in  nonuniform  macroscopic  3-D  domains.  Ap¬ 
plicability  to  nonunifonn  domains  is  a  clear  limitation  of 
existing  upscaled  mass  transfer  models  [e.g.,  Falta  et  al., 
2005a;  Wood  et  al.,  2005;  Basu  et  al.,  2006]. 

[10]  GTP  mass  ratios  quantify  the  distribution  of  DNAPL 
between  ganglia  and  pool  regions  according  to 

GTp  =  Y,Pnsn<t>AxAyAzds„  <  =  E  SrfSn  <  C**  ^ 

E  PnSnfpAxAyAzJs,,  >  i”ax  E  sn^sn  >  Cax 

Here  pooled  regions  are  defined  as  source  zone  regions 
(numerical  grid  blocks)  with  a  saturation  greater  than  the 


maximum  residual  organic  saturation  (s“ax),  which  corre¬ 
sponds  to  the  aqueous  phase  saturation  reached  upon 
complete  imbibition  from  the  residual  water  saturation  point 
on  the  drainage  curve  [Parker  and  Lenhard,  1987],  and 
ganglia  regions  are  defined  as  regions  with  DNAPL 
saturations  at  or  below  s“ax.  In  the  initial  stages  of  source 
evolution,  mass  discharge  is  controlled  by  the  dissolution  of 
high  surface  area  ganglia.  During  this  period,  GTP 
decreases  exponentially  with  time  (linearly  with  the  fraction 
of  mass  removed)  (results  not  shown).  At  later  time,  the 
persistence  of  DNAPL  pools  causes  the  source  zone  to 
transition  from  ganglia-dominated  to  pool-dominated.  Thus 
initial  source  zone  architecture  may  have  significant 
influence  on  mass  discharge  [Lemke  et  al.,  2004b;  Lemke 
and  Abriola,  2006]. 

[11]  Previous  research  on  upscaled  modeling  methods  has 
focused  on  the  simulation  of  DNAPL  source  zones  which 
are  characterized  in  this  work  as  end-members  on  a  contin¬ 
uum  of  potential  DNAPL  source  zone  distributions.  At  one 
end  of  the  continuum  of  source  zone  distributions  lies  the 
hypothetical  source  zone  with  a  uniform  DNAPL  distribu¬ 
tion  at  residual  saturation  (e.g.,  s„  =  0.01,  GTP  =  00,  used  by 
Zhu  and  Sykes  [2004]  to  simulate  mass  discharge).  Near  the 
opposite  end  of  the  continuum  lies  an  ensemble  of  hypo¬ 
thetical  source  zones  with  nonunifonn  DNAPL  distributions 
dominated  by  high  saturations  [Parker  and  Park,  2004; 
Park  and  Parker,  2005].  The  latter  distributions,  developed 
using  a  novel  percolation  theory-based  infiltration  model, 
resulted  in  relatively  high  entrapped  organic  phase  satura¬ 
tions  (i.e.,  contaminated  cell  average  sn  =  0.158;  GTP  <  1.0) 
in  the  upper  (U)  regions  of  the  domain  and  extensive 
pooling  above  the  lower  (L)  impermeable  boundary  (con¬ 
taminated  cell  average  sn  =  0.233;  GTP  <C  1.0),  which, 
when  combined  to  form  a  single  NAPL  distribution  (M) 
resulted  in  a  relatively  low  GTP  [Parker  and  Park,  2004; 
Park  and  Parker,  2005].  While  consideration  of  these 
scenarios  is  important,  field  conditions  will  likely  fall 
somewhere  between  the  ganglia-dominated  distribution  of 
Zhu  and  Sykes  [2004]  and  the  mixed  pool-dominated  source 
zone  distribution  of  Parker  and  Park  [2004]  [e.g.,  Poulsen 
and  Kueper,  1992;  Kueper  et  al.,  1993;  Broholm  et  al., 
1999], 

[12]  DNAPL  source  zones  employed  in  this  work  are 
composed  of  an  ensemble  of  16  3-D  DNAPL  saturation 
distributions  selected  from  a  recent  modeling  study  that 
simulated  hysteretic  infiltration  and  entrapment  of  an  im¬ 
miscible  tetrachloroethene  (PCE)-NAPL  in  nonunifonn 
permeability  fields  using  the  University  of  Texas  Chemical 
Flooding  Simulator  (UTCHEM  9.0)  [Center  for  Petroleum 
and  Geosystems  Engineering,  2000]  (see  Christ  et  al. 
[2005]  for  simulation  details).  Hydrogeologic  properties 
were  based  on  data  obtained  from  a  contaminated  aquifer 
in  Oscoda,  Michigan  (Bachman  Road  site),  composed  of 
relatively  homogeneous  glacial  outwash  sands  [Lemke  et 
al.,  2004a].  Spill  characteristics  were  selected  to  be  consis¬ 
tent  with  a  slow,  sustained  release  of  contaminant.  Simula¬ 
tion  parameters  are  summarized  in  Table  1.  Shown  in 
Figure  la  is  an  example  of  an  initial  DNAPL  saturation 
distribution.  Here,  the  simulation  domain  has  been  sliced 
through  the  center  cross  section  along  the  direction  of  flow  to 
show  the  PCE-NAPL  saturations  in  a  representative  x-y  plane. 
This  simulated  PCE-NAPL  distribution  is  consistent  with 
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Table  1.  Numerical  Simulation  Input  Parameters 


Parameter 

Water 

PCE 

Fluid  Properties 

Density  pj*  g  cirT 

0.999 

1.625 

Dynamic  viscosity,3  cP 

1.121 

0.89 

Compressibility,3  Pa-1 

4.4  x  10~‘° 

0.0 

Aqueous  diffusivity,b  cm2  s_I 

- 

8.6  x  10~6 

Aqueous  solubility,0  g  1 

- 

0.150 

Initial  saturation 

1.0 

0.0 

Spill  Scenario 

Spill  volume,  L 

128 

Spill  duration,  days 

400 

Release  rate,  L  m-2  d-1 

0.1725 

Pc-sa-krc  Model  Parameters 3 

Reference  air  entry  pressure,  kPa 

2.809 

Pore  size  index 

2.0773 

Interfacial  tension 

Air/water,  dyn  cm- 

72.75 

PCE/water,  dyn  cm-1 

47.8 

Irreducible  water  saturation 

0.080 

Max  residual  organic  saturation,  snrmax 

0.151 

Reference  permeability,  pm2 

19.7 

Variogram  Parameters 

Horizontal 

Vertical 

Matrix  Properties 3 

Nugget 

0.333 

0.333 

Range,  m 

7.0 

1.07 

Integral  scale,  m 

2.33 

0.36 

Mean  hydraulic  conductivity  K,a  m  d_  1 

16.8 

Anisotropy  ratio  kjkh 

0.5 

Lognormal  transformed  K  variance,3  cr2ln (K) 

0.29 

Applied  hydraulic  gradient,  m  m-1 

0.01 

Longitudinal  dispersivity  04, 3  m 

0.30 

Horizontal  transverse  dispersivity  ldhtT  m 

0.10 

Vertical  transverse  dispersivity  ujhv?  m 

0.0075 

Median  grain  size  d50,3  /mi 

295 

Uniformity  index  £/,-  3 

1.86 

Uniform  porosity  <pa 

0.36 

Ax  (m)  (Nx  =  26) 

0.3048 

Ay  (m)  (Ny  =  26) 

0.3048 

Az  (m)  (TV.  =  128) 

0.0726 

''Lemkc  et  al.  [2004a], 
bDekker  and  Abriola  [2000b], 
cHorvath  [1982], 
dUSEPA  [1986], 


those  generated  in  previous  modeling  studies  [e.g.,  Kueper 
and  Gerhard,  1995;  Dekker  and  Abriola,  2000a;  Lemke  et  al. , 
2004a]  and  observed  in  the  field  [e.g.,  Poulsen  and  Kueper, 
1992 \  Kueper  et  al.,  1 993 ;  Broholm  et  al. ,  1999].  In  contrast  to 
the  scenarios  used  by  Zhu  and  Sykes  [2004]  and  Park  and 
Parker  [2005],  the  DNAPL  spill  scenarios  simulated  herein 
resulted  in  a  mixture  of  low-saturation  ganglia  (0.0001  <5„< 
0.15)  and  high-saturation  pooled  (0.15  <sn<0  .7)  regions, 
without  accumulation  of  DNAPL  above  the  lower  imperme¬ 
able  boundary.  Hence  the  ensemble  of  saturation  distributions 
used  herein  is  distributed  between  two  end-members  repre¬ 
sented  by  the  existing  upscaled  modeling  studies  of  Zhu  and 
Sykes  [2004]  and  Parker  and  Park  [2004],  and  provides 
realistic  scenarios  along  the  continuum  of  potential  DNAPL 
source  zone  architectures. 

[13]  To  simulate  DNAPL  dissolution  and  mass  recovery 
from  the  selected  saturation  distribution  scenarios,  the  com¬ 
ponent  mass  balance  equation  was  solved  using  an  existing 
three-dimensional  (3-D),  multiphase  compositional  simula¬ 


tor;  the  modular  three-dimensional  transport  simulator 
(MT3D)  [Zheng  and  Wang,  1999],  modified  by  Parker  and 
Park  [2004].  Changes  in  bulk  aqueous  phase  flow  due  to 
DNAPL  dissolution  were  computed  using  MODFLOW 
(VERSION  1.10)  [McDonald  and  Harbaugh,  1988],  modi¬ 
fied  to  update  aqueous  phase  mobility  at  each  time  step 
[Parker  and  Park,  2004],  Interphase  mass  transfer  between 
the  aqueous  and  solid  phases  (i.e.,  sorption)  was  assumed  to 
be  negligible.  Mass  transfer  between  the  organic  liquid  and 
aqueous  phases  was  modeled  using  (2),  with  Kapi  obtained 
from  the  Powers  et  al.  [1994]  correlation: 


Ns.  =  4.13(ARe)0-598 


/  >  \  0.673 

(  _  ]  jyO.369 

\dM  J 


0.518+0.114  20  +0.10(7, 


(4) 


where  dM  is  the  diameter  of  a  “medium-sized”  sand  grain 
(dM  =  0.05  cm)  according  to  the  ASTM  particle  size  classi- 


S„ 

□.03 

0  025 

0.02 

0.015 

-0.01 

-0.005 


Figure  1.  (a)  Example  of  initial  DNAPL  saturation 
distribution  (GTP  =  23.0)  sliced  through  center  cross 
section  to  illustrate  nonuniform  DNAPL  saturations  in  x-y 
plane  and  (b)  flux-weighted  PCE  concentration  curves  (light 
gray  lines)  across  the  down-gradient  boundary  for  ensemble 
of  16  numerical  simulations.  Note  that  the  red  surface 
represents  the  sn  =  0.000 1  isosurface  and  that  the  range  of  ,v„ 
was  selected  to  differentiate  nonuniformities  in  low  s„ 
values.  Cells  with  s„  >  0.03  are  shaded  black. 
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fication,  Ut  is  the  uniformity  index  (c/60/<af i0),  s°n  is  the  DNAPL 
saturation  at  time  t  =  0,  and  all  other  parameters  are  as  defined 
previously.  The  modified  version  of  MT3D  was  used  to 
simulate  mass  recovery  until  99.99%  of  the  DNAPL  mass 
was  removed,  which  resulted  in  simulated  source  longevities 
between  5,000  and  ^30,000  days. 

[14]  Similar  correlations  have  been  used  extensively  in 
the  literature  [Mayer  and  Miller,  1996;  Powers  et  al,  1998; 
Zhu  and  Sykes,  2000;  Brusseau  et  al.,  2002]  to  describe  the 
local  (i.e.,  grid)  scale  mass  transfer  coefficient  as  a  function 
of  chemical  properties  and  system  parameters.  The  correla¬ 
tion  of  Powers  et  al.  [1994]  was  derived  using  effluent 
concentration  data  from  residual  NAPL  entrapped  in  one¬ 
dimensional  (1-D)  columns.  The  conditions  under  which  the 
Powers  correlation  was  developed  are  similar  to  those 
employed  in  the  majority  of  the  domains  simulated: 
0.0008  <  sn  <  0.181,  0.015  <  NRe  <  0.23,  and  1.19  <  Ut  < 
3.33.  Less  than  one  percent  of  the  DNAPL-contaminated 
cells  had  initial  saturations  exceeding  0.181.  Simulations 
conducted  with  alternative  correlations  suggest  that  the 
selection  of  a  mass  transfer  correlation  becomes  important 
only  at  very  late  time  (>95%  removal),  when  the  saturation 
distribution  is  typically  dominated  by  persistent  DNAPL 
pools  (results  not  shown).  Use  of  a  mass  transfer  correlation 
based  upon  ganglia  dissolution  in  this  phase  of  source 
evolution  (i.e.,  when  pools  dominate)  results  in  an  over 
prediction  of  down-gradient  concentrations.  While  this 
result  may  be  considered  conservative  for  the  prediction 
of  down-gradient  concentrations,  it  may  yield  an  under 
prediction  of  source  longevity.  This  is  discussed  further 
below  (see  section  3.4). 

2.3.  Analytical  Simulations 

[15]  Simplification  of  the  governing  equations  described 
in  section  2.1  is  required  to  derive  solutions  amenable  to 
analytical  techniques  [e.g.,  Sale  and  McWhorter,  2001;  Zhu 
and  Sykes,  2004;  Parker  and  Park,  2004].  Common  sim¬ 
plifying  assumptions  employed  in  upscaled  mass  transfer 
models  include  (1)  approximation  of  the  transient,  nonuni¬ 
form,  3-D  flow  field  by  a  steady,  uniform,  1-D  flow  field 
subject  to  an  equivalent  average  hydraulic  flux  ( q ),  (2)  neg¬ 
ligible  transverse  dispersivity  (ojt)  and  molecular  diffusion 
( Dm ),  (3)  quasi-steady  state  dissolution  kinetics,  and 
(4)  mass  transfer  (dissolution)  can  be  described  using  an 
upscaled  mass  transfer  coefficient  (ned)  that  is  a  function  of 
domain-averaged  parameters  [Dekker,  1996;  Parker  and 
Park,  2004;  Park  and  Parker,  2005].  Application  of 
assumptions  1  and  2  requires  that  the  down-gradient 
control  plane  is  sufficiently  large  such  that  all  of  the 
contaminant  mass  eluting  from  the  source  zone  crosses 
the  down-gradient  boundary.  In  practice,  the  down-gradient 
control  plane  must  increase  as  the  distance  separating  the 
source  zone  and  control  plane  increases  to  account  for 
mixing  in  directions  perpendicular  to  flow.  Assumption  3 
is  typically  justified  by  noting  that  the  timescale  for  changes 
in  DNAPL  saturations  (surrogate  for  interfacial  area)  during 
dissolution  is  relatively  large  when  compared  to  the  hy¬ 
draulic  residence  time  of  the  source  zone  [Parker  and  Park, 
2004], 

[16]  On  the  basis  of  the  first  assumption  (1,  uniform 
1-D  q),  the  bulk  fluid  phase  flow  equation  can  be  neglected. 


Application  of  assumptions  2-4  simplifies  the  component 
mass  balance  equation  (1)  to  yield 

(5) 

where  C(x)  is  the  flux-weighted  aqueous  phase  concentra¬ 
tion  of  a  selected  component  at  a  distance  x  from  the  source, 
where  x  is  typically  assumed  to  correspond  to  the  down- 
gradient  boundary  (x  =  L)  of  the  1-D  domain,  u,y  is  the 
longitudinal  dispersivity,  and  C’q  is  the  equilibrium  aqueous 
solubility  of  the  selected  component.  Equation  (5)  is  a 
second-order  ordinary  differential  equation  that  is  readily 
solved  by  enforcing  a  type  I  upstream  boundary  condition 
(C(0)  =  0),  a  type  II  boundary  condition  jylC(L)/dx  =  0)  at 
x  =  NxAx  =  L,  and  an  initial  condition  (C(x)  =  0)  [Parker 
and  Park,  2004]: 


C{L) 

Ceq 


=  1  —  exp 


L 

2  <j>u)L 


1-1 


Aneff(j)UL 

q 


which  is  readily  simplified  to 


C(L) 

Ceq 


(6a) 


(6b) 


for  uiJL  <0.1.  The  restriction  on  u>JL  is  reasonable,  given 
that  values  of  longitudinal  dispersivity  (04)  are  relatively 
small  when  compared  to  typical  distances  separating  a 
source  zone  from  the  selected  down  gradient  boundary  ( L ) 
(e.g.,  ujl  =  0.1  m  and  L  =  10  m). 

[17]  Application  of  (6a)  and  (6b)  requires  that  the 
upscaled  mass  transfer  coefficient  (k,,//)  be  quantified. 
Previous  studies  suggest  nef  may  be  evaluated  by  assuming 
that  it  conforms  to  a  field-scale  mass  transfer  correlation 
analogous  to  those  presented  for  local-scale  models  (e.g., 
(4)).  Best  fit  mass  transfer  correlation  parameters  have  been 
obtained  by  fitting  (6a)  and  (6b)  to  numerically  generated 
flux-weighted  concentration  curves  using  nonlinear  least 
squares  regression  [e.g.,  Parker  and  Park,  2004].  This 
fitting  process  yields  a  system-optimized,  upscaled  mass 
transfer  correlation,  which  approximates  the  dissolution 
behavior  of  a  particular  (but  unknown)  nonuniform  distri¬ 
bution  of  DNAPL  within  the  source  zone.  Application  of 
(6a)  and  (6b)  with  the  best  fit  upscaled  mass  transfer 
correlation  is  motivated  by  the  prospect  of  predicting  source 
zone  mass  discharge  at  other  sites  with  unknown  DNAPL 
saturation  distributions,  but  similar  source  zone  character¬ 
istics  [e.g.,  Parker  and  Park,  2004;  Falta  et  al.,  2005a;  Park 
and  Parker,  2005]. 

[is]  Source  longevity  predictions  may  also  be  made  using 
an  upscaled  mass  transfer  model  by  enforcing  a  mass 
balance  condition: 


—  =  - qAbC(L )  (7) 

Here  the  reduction  in  the  amount  of  DNAPL  mass  within 
the  source  zone  is  related  to  the  amount  of  mass  exiting  the 
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Table  2.  Summary  of  Simplified  Upscaled  Mass  Transfer  Models 


Reference 

Assumptions 

Regression  Equation 

Keff  Correlation 

Dekker  [1996] 

1,  2a 

(ks  dC_ 
Wa  Qt 

|(<M“  f)+?S  =  v/(ce*-c) 

Z  (  r. aVS  \  P 

zpen  1  Sa  \ 

NzAz  J 

solved 

using  finite  element  method 

K°  (l)  (m(/=0))  7  “  1  m  all  simulations 

Parker  and  Park  [2004] 

1,  2,  3,  4a 

note 

«  ^  for  'vAy-  <  1 

Zhu  and  Sykes  [2004]b 

CaM(t) 

.  C(L)  _  M(t) 

C«7  —  M(t= 0) 

Keffno\.  specified  in  this  model 

NL1.c(u_  f  m  Y 

NL0-  C{L)  —  /  M{t)  \  ^ 

INLZ.  ceq  -  K0yM(t= 0) J 

Falta  et  al.  [2005 a]c 

CaM(t) 

C(L)  _  (  M(t)  Y 

Co  \M{t= 0)  J 

Kejj  not  specified  in  this  model 

Assumption  described  in  text. 

bZhu  and  Sykes  [2004]  present  three  simplified  models:  L,  linear  model;  NL1,  nonlinear  model  1;  and  NL2,  nonlinear  model  2. 
CC„  reflects  the  flow-averaged  contaminant  concentration  at  the  source  zone  control  plane. 


domain  through  the  down  gradient  boundary  ( Ab )  in  (7). 
Substitution  of  expressions  for  C(L)  derived  from  the 
upscaled  models  into  (7)  provides  a  simple  method  for 
computing  the  time  required  to  achieve  a  specified  mass 
reduction  based  on  remediation  objectives. 


assuming  a  basic  relationship  between  the  flux-weighted 
effluent  concentration  and  mass  depletion  [Zhu  and  Sykes, 
2004;  Falta  et  al.,  2005a].  Regardless  of  model  complexity, 
however,  all  upscaled  mass  transfer  correlations  employ  the 
same  functional  form: 


3.  Results 

3.1.  Baseline  Numerical  Simulations 


(8) 


[19]  As  described  in  section  2.2,  an  ensemble  of  16 
numerical  simulations  of  mass  discharge  from  a  nonuniform 
PCE  source  zone  were  developed  and  used  as  a  baseline  for 
comparison  with  predictions  obtained  from  four  existing 
upscaled  mass  transfer  models  [Parker  and  Park,  2004;  Zhu 
and  Sykes,  2004].  Flux-weighted  aqueous  phase  concen¬ 
trations  across  the  down-gradient  domain  boundary  for  all 
16  realizations  are  plotted  in  Figure  lb.  The  assumption  that 
PCE-NAPL  infiltration  occurs  in  the  absence  of  dissolution 
results  in  an  initial  flux-weighted  concentration  across  the 
down-gradient  domain  boundary  equal  to  zero.  This  flux- 
weighted  concentration  rises  sharply  at  very  low  percent 
mass  removal  as  dissolution  and  transport  of  contaminant 
begins.  Down-gradient  concentrations  reach  “pseudo¬ 
steady  state”  values  (^40  mg  L“')  after  approximately  5% 
mass  removal.  Although  saturation  distributions  change  due 
to  DNAPL  dissolution,  the  resulting  change  in  flux-weighted 
concentrations  is  relatively  small  until  the  source  zone  begins 
to  transition  from  a  ganglia-dominated  to  a  pool-dominated 
DNAPL  distribution.  This  transition  region  occurs  between 
75  and  95%  mass  removal  (1,000-2,000  days)  in  these 
simulations,  and  results  in  a  rapid  reduction  of  contaminant 
concentration.  Concentrations  approach  zero  only  after  the 
vast  majority  (>95%)  of  the  contaminant  mass  has  been 
removed  from  the  persistent  high  saturation  pools. 

3.2.  Upscaled  Model  Assessment 

[20]  The  general  method  used  in  upscaled  mass  transfer 
modeling  (section  2.3)  has  been  employed  in  the  literature 
along  with  various  simplifying  assumptions  to  derive  sev¬ 
eral  upscaled  mass  transfer  models  for  specific  source  zone 
scenarios.  These  models  are  summarized  in  Table  2.  In¬ 
spection  of  Table  2  shows  that  the  existing  upscaled  models 
range  in  complexity  from  a  nonsteady,  numerically  based 
model  [Dekker,  1996]  to  highly  simplified  models  derived 


where  k0  and  /3  are  fitting  parameters  and  Ma  is  the 
DNAPL  mass  at  t  =  0.  Here,  k'0  is  a  function  of  system 
parameters  (q  and  L)  and  is  used  to  account  for  mixing 
and  dilution  as  water  flows  through  the  nonuniform 
DNAPL  saturation  distribution  (i.e.,  mass  discharge).  Park 
and  Parker  [2005]  recently  presented  a  methodology  for 
computing  k0  based  on  aquifer  and  source  zone  properties. 
In  contrast,  (3  is  a  fit  parameter  used  as  a  surrogate  for  the 
saturation  architecture  that  controls  the  effective  rate  at 
which  mass  transfer  occurs  during  depletion  of  DNAPL 
mass. 

[21]  Sensitivity  of  these  two  parameters  to  source  zone 
conditions  can  be  explored  by  comparing  values  derived  for 
different  source  zone  scenarios  (see  Table  3).  To  develop 
row  1  of  Table  3,  the  upscaled  model  (6a)  and  (6b)  was  fit  in 
conjunction  with  the  upscaled  mass  transfer  correlation  (8) 
to  each  realization  of  the  ensemble  of  numerically  generated 
mass  recovery  curves  (see  Figure  lb).  Mean  best  fit 
upscaled  correlation  parameter  values  and  their  cor¬ 
responding  95%  confidence  intervals  are  presented  in  the 
first  row  of  Table  3.  The  flux-weighted  concentration 
prediction  obtained  from  the  mean  of  the  best  fit  values  in 
the  upscaled  mass  transfer  model  is  compared  to  the 
ensemble  of  numerical  simulations  in  Figure  2.  For  com¬ 
parison  purposes,  Table  3  also  presents  mass  transfer 
parameters  reported  for  the  correlations  presented  in 
Table  2.  Recall  that  the  correlations  shown  in  Table  2  were 
developed  from  source  zone  scenarios  which  differ  from 
those  based  on  the  Bachman  Road  site.  Hence  these 
correlations  are  not  expected  to  result  in  accurate  flux- 
weighted  concentration  predictions  for  the  Bachman  Road 
site.  Predictions  of  flux-weighted  concentrations  shown  in 
Figure  2  illustrate  the  variability  that  can  be  expected  when 
an  existing  correlation  is  inappropriately  applied  to  a  site 


6  of  13 


W11420 


CHRIST  ET  AL.:  ESTIMATING  DISCHARGE  FROM  DNAPL  SOURCE  ZONES 


W11420 


Table  3.  Comparison  of  Best  Fit  and  Published  Correlation  Parameters  for  Ke^-in  Upscaled  Mass  Transfer  Models  Shown  in  Table  2a 


Reference 

System  Description 

<  cT1 

0 

This  work 

(equations  (6)  and  (8)) 

ensemble  of  3-D  simulations  of  nonuniform  source  zones 
based  on  Bachman  Road  site 

8.2  x  1(T3  (±0.06) 

0.85  (±0.07) 

Dekker  [1996] 

2-D  simulation  of  SEAR  in  nonuniform  DNAPL 

source  zones  based  on  Borden  (B)  and  Jussel  (J)  aquifers’3 

B:  0.14  J,  0.09 

B:  0.81  J,  1.5 

Parker  and  Park  [2004] 

3-D  simulation  of  DNAPL  dissolution 
under  natural  gradient  conditions  in  upper  region  (U), 
lower  region  (L),  and  upper  and  lower  region  (M) 
of  a  single  nonuniform  DNAPL  source  zone 

U:  4.9  x  10~4  L,  2.4  x 
10~s  M,  4.9  x  10~4 

U:  1.10  L,  0.40  M,  1.40 

Zhu  and  Sykes  [2004] 

1-D  column  containing  a  residual  saturation  of 

DNAPL  using  the  linear  (L),  nonlinear  1  (NL1), 
and  nonlinear  2  (NL2)  models0 

L:  NA  NL1,  NA  NL2,  0.80 

L:  NA  NL1,  1.10  NL2,  1.31 

Falta  et  al.  [2005a] 

3-D  simulation  of  NAPL  dissolution 
under  natural  gradient  conditions  with  organic  saturation 
positively  or  negatively  correlated  with  permeability 

NA 

0. 5-2.0 

aNA,  not  applicable. 

bThe  Borden  and  Jussel  aquifers  differ  in  their  degree  of  spatial  heterogeneity:  a2ln(k)  =  0.24  for  Borden  and  ir2ln(A')  =  1 .0  for  Jussel  [see  Dekker,  1 996; 
Dekker  and  Abriola,  2000a], 

cZhu  and  Sykes  [2004]  state  that  this  comparison  is  for  illustrative  purposes  only  and  was  not  a  true  test  of  the  models. 


with  characteristics  that  are  not  consistent  with  those  to 
which  the  correlation  is  calibrated.  Note,  the  Dekker  [1996] 
model  is  excluded  from  Figure  2  because  this  model 
employs  numerical  methods,  which  diminishes  a  key  benefit 
of  upscaled  modeling  approaches;  ease  of  use.  The  Falta  et 
al.  [2005a]  model  has  also  been  excluded  from  Figure  2 
because  it  is  essentially  the  same  model  as  the  Zhu  and 
Sykes  [2004]  NL2  model. 

[22]  Inspection  of  the  parameter  values  reported  in  Table  3 
reveals  differences  among  model  parameters  arise  from 
differences  in  the  source  zone  conditions  used  to  obtain 
the  best  fit  correlation  parameters.  For  example,  Parker  and 
Park  [2004]  simulated  the  release  of  229  L  of  trichloroe- 
thene  (TCE)  with  a  release  rate  of  ~25,000  L  m  2  d”1  from 
a  single  cell  into  a  randomly  generated  permeability  field 
with  a  mean  saturated  hydraulic  conductivity  ( K )  of  10  m  d_  1 
and  a  ln(AT)  variance  (cr2ln (K))  of  1.  This  relatively 
high  release  rate,  combined  with  a  lack  of  capillary  entrap¬ 
ment  in  the  percolation  model,  led  to  high  TCE  saturations 
(GTP  «  0.414  assuming  “fingers”  and  “lenses”  values 
presented  in  Table  1  of  Parker  and  Park  [2004]  correspond 
to  ganglia  and  pools,  respectively)  distributed  nonuniformly 
throughout  a  relatively  small  region  of  the  domain  (~10% 
of  the  10  x  10  m  footprint)  and  resulted  in  dilution  of  the 
flux-weighted  concentration  across  the  down-gradient 
boundary  (C(Z)).  Fitting  this  diluted,  flux-weighted  concen¬ 
tration  resulted  in  n'a  values  that  are  more  than  an  order 
of  magnitude  less  than  the  parameters  obtained  with  the 
other  correlations  (see  Table  3),  regardless  of  which  region 
of  the  domain  was  considered  (upper  (U),  lower  (L),  or 
mixed  (M),  see  Table  2).  As  shown  in  Figure  2  for  the 
conditions  considered  here,  a  low  k0  value  translates  to  a 
relatively  low  quasi-steady  state  contaminant  concentration 
prediction.  In  contrast,  Zhu  and  Sykes  [2004]  considered  a 
hypothetical  2-D  (x-z)  domain  that  contained  a  fully  pene¬ 
trating  (z  direction),  residually  saturated  ( sn  =  0.01)  region 
with  a  GTP  =  00.  This  relatively  simple  source  zone 
scenario  leads  to  higher  best  fit  k0  values,  which  result  in 
near-saturation  quasi-steady  state  contaminant  concentration 
predictions  (Figure  2).  Variability  among  the  fitted  values  of 
j3  (mass  depletion  parameter)  was  within  a  factor  of  three. 
Low  values  of  /3  (<1.0)  generally  correspond  to  saturation 


distributions  dominated  by  DNAPL  pools.  High  (3  values 
(>1.0)  correspond  to  DNAPL  saturation  distributions  dom¬ 
inated  by  DNAPL  ganglia. 

[23]  Relative  differences  between  upscaled  model  predic¬ 
tions  of  flux-weighted  concentrations  and  the  mean  flux- 
weighted  concentration  of  the  numerical  simulations  can  be 
quantified  using  objective  functions  such  as  the  root-mean- 
square  error  (RMSE).  RMSE  was  computed  at  5%  mass 
removal  increments  starting  at  the  mass  removal  level 
corresponding  to  the  establishment  of  quasi-steady  state 
conditions  (i.e.,  5%)  and  ending  at  99.99%  mass  removal. 
RMSE  values  for  the  upscaled  models  depicted  in  Figure  2 
are  listed  in  Table  4.  Use  of  an  upscaled  model  developed 
under  alternative  site  conditions  to  predict  flux-weighted 
concentrations  for  the  conditions  of  the  numerical 
simulations  considered  in  this  work  results  in  either  an 
overprediction  ( Zhu  and  Sykes  [2004]  models)  or  under¬ 
prediction  ( Parker  and  Park  [2004]  M)  of  the  flux-weighted 
concentration.  The  upscaled  model  fit  to  the  ensemble  of 
16  numerical  simulations  provides  the  lowest  RMSE.  It  is 
important  to  note,  however,  that  even  when  the  upscaled 


Figure  2.  Flux-weighted  PCE  concentration  as  a  function 
of  mass  reduction  for  16  numerical  simulations  (light  gray 
lines)  and  upscaled  mass  transfer  models  listed  in  Table  2. 
Model  input  parameters  are  given  in  Tables  1,  2,  and  3. 
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Table  4.  Root-Mean-Square  Error  of  Upscaled  Models  When 
Compared  to  the  Mean  of  the  Numerical  Simulations3 


Upscaled  Model 

RMSE,  mg  L_I 

This  work  (equations  (6)  and  (8)) 

1.8 

Parker  and  Park  [2004]  M 

27.3 

Zhu  and  Sykes  [2004]  L 

53.0 

Zhu  and  Sykes  [2004]  NL1 

50.1 

Zhu  and  Sykes  [2004]  NL2 

30.6 

“Note  that  RMSE  computed  at  mass  removal  increments  of  5%  starting  at 
the  5%  mass  removal  level  and  ending  at  99.99%  mass  removal. 


mass  transfer  parameters  are  fit  to  the  numerical  simulations 
(this  work),  flux-weighted  concentration  predictions  can 
overestimate  or  underestimate  the  flux-weighted  concentra¬ 
tion  for  any  individual  realization  (light  gray  lines  in 
Figure  2)  by  an  order  of  magnitude. 

[24]  An  alternative  metric  for  evaluating  the  perfonnance 
of  upscaled  mass  transfer  models  is  source  longevity  (equa¬ 
tion  7).  Source  longevity  predictions  using  upscaled  models 
are  compared  to  source  longevity  statistics  for  the  16 
numerical  simulations  in  Table  5.  Overestimation  of  the 
flux-weighted  concentrations  [e.g.,  Zhu  and  Sykes,  2004] 
resulted  in  relatively  short  cleanup  times  (reduced  source 
longevity),  whereas  underestimation  of  flux-weighted  con¬ 
centrations  (e.g.,  Parker  and  Park  [2004]  M  model)  resulted 
in  increased  source  longevity  (by  up  to  one  order  of 
magnitude).  The  time  required  to  achieve  higher  levels  of 
mass  removal  (e.g.,  99.99%)  has  been  shown  to  depend  on 
the  presence  or  absence  of  persistent,  discontinuous 
DNAPL  pools  [Lemke  et  al.,  2004b;  Lemke  and  Abriola, 
2006],  Given  that  the  mass  depletion  exponent,  (3,  is  derived 
by  fitting  the  entire  DNAPL  mass  depletion  curve,  it  is 
unlikely  that  upscaled  mass  transfer  models  can  predict 
flux-weighted  concentrations  at  large  times,  when  DNAPL 
ganglia  have  dissolved  and  discrete  DNAPL  pools  remain. 
Hence,  as  the  level  of  mass  removal  increases  from  90%  to 
99.99%,  the  disparity  between  upscaled  source  longevity 
predictions  and  the  average  source  longevity  prediction 
obtained  from  the  numerical  simulations  will  likely  grow 
larger. 

3.3.  Pseudosteady  State  Concentration 

[25]  The  analysis  presented  above  indicates  that  predic¬ 
tions  of  flux-weighted  concentrations  and  source  longevity 
using  upscaled  models  can  substantially  deviate  from  nu¬ 
merical  simulation  results  when  the  models  are  applied  to 


sites  with  NAPL  saturation  distributions  uncharacteristic  of 
those  that  were  used  to  calibrate  the  upscaled  correlation 
parameters.  Most  upscaled  model  predictions  shown  in 
Figure  2  result  in  higher  flux-weighted  concentrations  than 
those  obtained  for  the  average  numerical  simulation,  which 
may  be  considered  a  conservative  estimate  with  respect  to 
maximum  contaminant  level  (MCL)  compliance  criteria. 
However,  higher  flux-weighted  concentration  values  also 
result  in  the  underestimation  of  source  longevity  or  cleanup 
time  (see  Table  5).  It  may  therefore  be  advantageous  to 
modify  the  upscaled  mass  transfer  correlation  (8)  to  obtain 
more  accurate  estimates  of  the  initial  flux-weighted  con¬ 
taminant  concentrations. 

[26]  Two  possible  approaches  for  improving  the  predic¬ 
tive  capabilities  of  the  upscaled  mass  transfer  correlation  are 
( 1 )  incorporation  of  additional  parameters  in  the  correlation 
to  account  for  the  effects  of  flow  bypassing  on  the  flux- 
weighted  concentrations  and  (2)  modification  of  an  existing 
parameter  in  the  correlation  (e.g.,  na)  to  more  accurately 
represent  flux-weighted  concentrations.  Dekker  [1996] 
employed  the  first  type  of  modification  by  incorporating 
the  depth  of  the  DNAPL  saturation  distribution  (depth  of 
penetration,  zpen)  normalized  by  the  depth  to  the  lower 
impermeable  boundary  ( NzAz )  into  the  upscaled  mass 
transfer  correlation  to  account  for  dilution  at  the  down- 
gradient  boundary  (see  Table  2).  An  analogous  approach 
could  be  used  in  3-D  to  modify  the  down  gradient  boundary 
area  to  be  consistent  with  the  horizontal  spread  of  the 
saturation  distribution,  e.g., 


Ac 

Keff  ~  Ko  ^ 


(  M(t)  y 

\M(‘  =  0  )J 


(9) 


where  Ac  is  the  contaminated  area  of  the  down  gradient 
boundary  area  ( Ah ).  Quantification  of  Ac  requires  delinea¬ 
tion  of  the  extent  of  the  dissolved-phase  contaminant  plume 
at  a  selected  concentration  (e.g.,  MCL).  Precise  delineation 
of  the  contaminant  concentration  profile  is  a  potentially 
difficult  task  in  the  field,  particularly  for  source  zones  with 
high  degrees  of  nonuniformity. 

[27]  Alternatively,  inspection  of  the  best  fit  upscaled 
simulation  shown  in  Figure  2  indicates  that  the  initial 
flux-weighted  concentration  (t  «  0,  M/M0  «  1)  is  approx¬ 
imately  equivalent  to  the  pseudo-steady  state  flux-weighted 
concentration  predicted  using  the  upscaled  model.  Assum¬ 
ing  that  the  initial  flux-weighted  concentration  (CQ)  could 
be  obtained  from  preliminary  site  characterization  data, 


Table  5.  Source  Longevity  Predicted  Using  the  Numerical  and  Upscaled  Mass  Transfer  Models3 


Numerical  Simulation 

Parker  and  Park  [2004]  M 

Zhu  and  Sykes  [2004] 

Upscaled  Model  (6) 

Modified  Upscaled  Model  (14) 

Minimum  2.2 

Mean  4.9 

Maximum  12.7 

130.2 

90%  Mass  Removal 

L,  0.9 

NL1,  1.9 

NL2,  2.6 

minimum  4.4 
mean  5.7 
maximum  8.8 

minimum  4.4 
mean  5.6 
maximum  7.2 

Minimum  13.9 

Mean  27.4 

Maximum  89.7 

1279 

99.99%  Mass  Removal 

L,  3.7 

NL1,  4.8 

NL2,  6.7 

minimum  8.3 
mean  20.8 
maximum  102.2 

minimum  7.6 
mean  14 
maximum  29.7 

“Source  longevity  is  defined  here  as  the  time  in  years  to  remove  90%  or  99.99%  of  the  initial  DNAPL  mass. 
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Percent  Mass  Removal 

Figure  3.  Flux-weighted  PCE  concentrations  as  a  function 
of  mass  reduction  for  16  numerical  simulations  (light  gray 
lines)  and  upscaled  mass  transfer  models  modified  as 
described  in  the  text  (equations  (11)  and  (12)).  Model 
parameters  are  given  in  Tables  1  and  3. 


equations  reported  in  Table  2  may  be  rearranged  to  com¬ 
pute,  rather  than  fit,  k0.  For  example,  the  Parker  and  Park 
[2004]  equation  can  be  rearranged  to  compute  nB  as  a 
function  of  system  properties  and  C0: 

K“  =  -zln(1'^)  (10) 

Flux-weighted  concentrations  at  the  down-gradient  bound¬ 
ary  may  then  be  determined  by  substituting  (10)  into  the 
regression  equation  used  to  obtain  best  fit  parameters  for  Keff 
((6a),  (6b),  and  (8)): 


values  are  consistent  with  the  RMSE  obtained  when  fitting 
the  upscaled  mass  transfer  model  (equations  6a,  6b  and  8)  to 
the  numerical  simulations,  but  require  only  one  measured 
(CQ)  and  one  best  fit  ((3)  parameter. 

3.4.  Source  Zone  Architecture 

[29]  A  number  of  researchers  have  recognized  the  need  to 
incorporate  site-specific  information  describing  the  source 
zone  architecture  into  upscaled  modeling  predictions 
[Dekker,  1996;  Parker  and  Park,  2004;  Falta  et  ai, 
2005a].  As  shown  in  the  previous  analysis,  this  information 
is  necessary  to  accurately  predict  the  rate  at  which  DNAPL 
mass  transfer  occurs  in  the  source  zone  and  thus  may 
influence  the  mass  depletion  parameter  (3).  In  recent 
numerical  modeling  studies,  Lemke  et  al.  [2004b]  and 
Lemke  and  Abriola  [2006]  showed  that  the  distribution  of 
DNAPL  between  ganglia  and  pool  regions,  defined  using 
the  GTP  mass  ratio  (see  section  2.2),  controlled  the  flux- 
weighted  concentration  emanating  from  the  source  zone.  It 
is  hypothesized  here  that  the  GTP  metric  is  a  better 
predictor  of  DNAPL  mass  dissolution  characteristics  in  a 
source  zone  than  more  traditional  metrics  (e.g.,  center  of 
mass,  average  DNAPL  saturation). 

[30]  A  relationship  between  the  mass  depletion  exponent 
(/?)  and  the  initial  GTP  saturation  distribution  was  devel¬ 
oped  from  an  analysis  of  the  1 6  numerical  simulations  used 
in  this  study.  Flux-weighted  concentration  curves  for  each 
simulation  were  fit  to  (11)  by  setting  Ca  equal  to  the 
numerically  simulated  pseudo-steady  state  flux-weighted 
concentration  and  then  fitting  /3  using  a  nonlinear  least 
squares  regression  routine.  Shown  in  Figure  4  is  a  plot  of  (3 
versus  the  GTP  for  each  of  the  numerical  simulations.  A 
statistically  significant  (p  value  =  0.0039)  negative  correla¬ 
tion  was  obtained  between  f3  and  GTP,  which  can  be 
expressed  through  a  power  law  model  (r1  =  0.6373)  in  the 
following  form: 


C(L) 


(11) 


Only  one  fitting  parameter  ( (3 )  is  included  in  (11).  The  best 
fit  parameter,  has  been  replaced  with  a  parameter  (Cc) 
that  may  be  measured  in  the  field.  Perfonning  this  same 
analysis  for  the  other  upscaled  mass  transfer  models  (Zhu 
and  Sykes  [2004],  Falta  et  al.  [2005a],  and  simplified  (i.e., 
KeffL!q  C  1)  Parker  and  Park  [2004])  yields 


C{L)  =  C0  f M\ 
C e<i  Ce«  \Mo  J 


(12) 


Note  that  the  dependence  of  equation  (12)  on  the  source 
zone  architecture  (/?),  despite  the  assumption  that  KeffLlq  <C 
1  (i.e.,  nonequilibrium),  stems  from  the  assumption  that  neff 
is  linearly  related  to  q/L. 

[28]  Improvements  in  the  predictions  obtained  using  (11) 
or  (12),  depending  on  the  model,  are  shown  in  Figure  3.  All 
upscaled  mass  transfer  models  now  yield  pseudo-steady 
state  flux-weighted  concentration  predictions  within 
67%  of  the  mean  of  the  numerical  simulations.  The 
corresponding  RMSE  for  each  of  the  modified  models 
was  reduced  by  a  factor  between  3  and  10.  These  RMSE 


j3=  1.5  ■  GTP-026  for  1.5  <  GTP  <  24.0  (13) 

Incorporation  of  (13)  into  (11)  and  (12)  results  in  an 
upscaled  mass  transfer  model  capable  of  simulating  both  the 


Figure  4.  Mass  depletion  parameter  (J3)  as  a  function  of 
ganglia-to-pool  (GTP)  mass  ratio  for  15  numerical 
simulations.  One  simulation  with  GTP  =  00  was  excluded 
from  the  figure  and  correlation. 
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pseudo-steady  state  flux-weighted  concentration  and  the 
transition  from  ganglia-  to  pool-dominated  dissolution: 


_  _  ( M 

C(L )  f  C0\l* 

Ceq  ^  Ce1 ) 


(14a) 


C{L) 

Ceq 


c^(m  y-5^' 

c *1  \M0) 


(14b) 


Average  contaminant  concentrations  as  a  flinction  of  two 
site-specific  parameters,  the  initial  flux-weighted  concen¬ 
tration  (CQ)  and  the  source  zone  GTP,  are  included  in  (14a) 
and  (14b).  This  formulation  does  not  require  fitting  to 
experimental  or  numerical  simulations  in  order  to  obtain 
upscaled  correlation  parameters.  Rather,  site-specific  field 
measurements  (i.e.,  C0  and  GTP)  are  incorporated  directly 
into  the  upscaled  model,  providing  a  more  versatile  method 
for  the  prediction  of  the  flux-weighted  concentration 
corresponding  to  a  given  level  of  mass  removal  (M/M0). 

[31]  Comparing  the  proposed  upscaled  model  (14a)  pre¬ 
diction  of  flux-weighted  concentration  to  the  mean  of  the 
numerical  simulations  resulted  in  a  RMSE  of  approximately 
2  mg  L  1 .  This  RMSE  is  almost  identical  to  the  RMSE  for  the 
original  upscaled  mass  transfer  model  (1.8  mg  L-1,  Table  4). 
The  mean  source  longevity  value  for  the  proposed  (14a)  and 
original  (6a)  and  (6b)  upscaled  mass  transfer  models  is  nearly 
identical  at  90%  mass  removal  (<2%  difference),  and  within 
50%  at  99.99%  mass  removal  (Table  5).  These  comparisons 
demonstrate  the  ability  of  the  proposed  upscaled  model  to 
reproduce  the  predicted  values  obtained  using  the  original 
best  fit  upscaled  model  without  the  need  to  fit  scenario- 
specific  correlation  parameters. 

[32]  To  examine  the  ability  of  the  model  to  be  used  in  a 
predictive  sense,  however,  the  proposed  model  (14a)  and 
(14b)  must  be  applied  to  source  zone  scenarios  other  than 
those  used  to  develop  the  correlation  (13).  For  this 
evaluation,  an  ensemble  (10  realizations)  of  published 
source  zone  scenarios  was  selected.  Eight  of  the  source 
zone  realizations  were  selected  from  a  recent  numerical 
modeling  study  that  investigated  the  relationship  between 
mass  discharge  and  source  zone  architecture  [Lemke  and 
Abriola,  2006].  These  realizations  were  developed  using 
sequential  indicator  simulation  (SIS)  and  assuming  there 
was  no  correlation  between  soil  (i.e.,  permeability)  and 
capillary  properties  (i.e.,  entry  pressure).  Source  zone 
saturation  distributions  generated  using  SIS  were  charac¬ 
terized  by  high  saturations  and  extensive  lateral  pooling. 
GTP  mass  ratios  ranged  from  0.9  to  65  and  mass  recovery 
was  simulated  assuming  a  surfactant  enhanced  aquifer 
remediation  scenario.  The  remaining  two  realizations  cor¬ 
respond  to  the  simulations  used  by  Zhu  and  Sykes  [2004] 
and  Parker  and  Park  [2004],  which  have  GTP  mass  ratios 
of  infinity  and  ^0.41,  respectively.  Source  zones  with 
GTP  values  outside  the  range  used  in  the  development 
of  the  correlation  were  selected  to  show  the  limitations  of 
the  model  at  low  or  high  GTP  values.  Flux-weighted 
concentration  predictions  made  using  the  proposed 


upscaled  model  (solid  line)  are  compared  to  the  numerical 
simulation  results  (open  symbol)  in  Figure  5.  GTP  mass 
ratios  and  RMSE  for  each  source  zone  simulation  are  also 
reported  in  Figure  5.  Initial  flux-weighted  concentrations 
(C0)  were  assumed  to  correspond  to  the  first  numerically 
simulated  concentration  value.  For  the  simulations  of 
Lemke  and  Abriola  [2006]  this  value  did  not  necessarily 
correspond  to  a  0%  mass  removal  level,  but  was  assumed 
to  be  a  reasonable  approximation  to  the  initial  flux- 
weighted  concentration  value. 

[33]  Analysis  of  Figure  5  demonstrates  that  the  proposed 
upscaled  model  predicts  flux-weighted  concentrations 
with  RMSE  values  that  range  from  3  to  20  mg  L_I. 
These  RMSE’s  are  generally  three  to  five  times  less  than 
the  RMSE  values  obtained  when  the  previously  published 
models  were  used  for  prediction  (Table  4).  The  upscaled 
model  captured  the  initial  pseudo-steady  state  flux-weighted 
concentration  and  generally  captured  the  transition  from  a 
ganglia-dominated  to  a  pool-dominated  source  zone,  but  in 
some  cases  failed  to  accurately  predict  the  flux -weighted 
concentration  eluting  from  the  source  zone  at  high  levels 
(>80%)  of  mass  removal  (e.g.,  Figure  5c).  As  discussed 
previously,  this  limitation  stems  from  the  assumption  that  the 
mass  depletion  exponent  is  a  constant  and  that  the  correla¬ 
tion  is  based  upon  scenarios  characterized  by  relatively  high 
GTP  mass  ratios.  Substituting  (13),  which  is  a  function  of  the 
initial  GTP  mass  ratio,  for  this  mass  depletion  exponent 
improves  the  ability  of  the  upscaled  model  to  simulate  the 
transition  from  a  ganglia-  to  a  pool-dominated  source  zone, 
however,  it  does  not  account  for  the  eventual  transition  to  a 
source  zone  dominated  by  discreet  DNAPL  pools.  The 
observation  that  this  discrepancy  generally  occurs  for  source 
zones  with  low  initial  GTP  mass  ratios  (Figures  5a-5c)  is  a 
further  indication  of  the  importance  of  DNAPL  pools  in 
determining  the  contaminant  flux  eluting  from  a  nonuniform 
source  zone. 

[34]  A  second  limitation  of  this  model  is  the  need  to 
estimate  the  initial  DNAPL  source  zone  GTP  mass  ratio. 
Currently,  there  is  a  lack  of  methods  to  accurately  quantify 
DNAPL  saturation  distributions  in  the  field.  Although 
partitioning  interwell  tracer  tests  (PITT)  have  been  used 
to  estimate  average  NAPL  saturations  in  source  zones  over 
the  integrated  path  length  [e.g.,  Amiable  et  al.,  1998; 
Ramsburg  et  ah,  2005],  the  results  are  spatially  averaged 
(e.g.,  1  m)  and  thus  do  not  provide  the  resolution  required 
to  distinguish  between  ganglia-dominated  and  pool- 
dominated  regions.  Resolution  in  NAPL  saturation  esti¬ 
mates  resulting  from  PITTs  may  be  increased  by  locating 
multiple  observation  points  along  a  flow  path  [Jin  et  al., 
1995;  Annable  et  ah,  1998;  Rao  et  ah,  2000].  Pure  et  ah 
[2006]  suggest  that  higher-order  temporal  moments 
employed  for  the  interpretation  of  PITT  data  [Jawitz  et 
ah,  2003]  may  hold  promise  for  identifying  source  zone 
metrics  which  link  mass  discharge  to  NAPL  architecture. 
In  addition,  it  may  be  possible  to  obtain  higher  resolution 
estimates  of  the  spatial  distribution  of  DNAPL  mass 
through  the  use  of  geophysical  techniques,  which  have 
been  used  to  directly  estimate  in  situ  DNAPL  saturations 
[e.g.,  Chambers  et  ah,  2004],  The  analysis  presented 
herein  clearly  demonstrates  the  importance  of  incorpora¬ 
tion  of  site-specific  parameters  into  upscaled  mass  transfer 
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Figure  5.  Comparisons  of  numerical  simulation  (open  symbols)  and  proposed  upscaled  mass  transfer 
model  (14a)  predictions  (solid  lines)  of  flux- weighted  PCE  concentrations  as  a  function  of  mass  removal 
for  realizations  selected  from  Lemke  and  Abriola  [2006]  (Figures  5a-5h)  and  Zhu  and  Sykes  [2004] 
(Figure  5i)  source  zone  scenario  (GTP  =  oo),  and  Parker  and  Park  [2004]  source  zone  scenario  (GTP  = 
0.414).  GTP  mass  ratios  and  RMSE  are  (a)  1.85  and  15.7  mg  L_1,  (b)  0.91  and  9.0  mg  L_1,  (c)  1.78  and 
10.7  mg  L_1,  (d)  6.37  and  10.0  mg  L_1,  (e)  1.75  and  6.4  mg  L-1,  (f)  3.18  and  7.7  mg  L-1,  (g)  65.6  and 
19.5  mg  L_1,  (h)  30.6  and  10.0  mg  L_1,  and  (i)  oo  and  3.05  mg  L  1  for  Zhu  and  Sykes  [2004]  and  0.41 
and  3.4  mg  L-1  for  Parker  and  Park  [2004],  respectively. 


models  and  highlights  the  need  for  innovations  in  in  situ 
characterization  techniques  as  an  important  area  for  future 
research. 

4.  Summary  and  Conclusions 

[35]  There  is  a  growing  need  to  develop  relatively 
simple  modeling  tools  to  aid  site  managers  and  regulators 
in  the  monitoring  and  treatment  of  DNAPL  contaminated 
source  zones.  Site-specific  numerical  simulations  are 
among  the  most  effective  methods  for  evaluating  alterna¬ 
tive  management  strategies;  however,  such  models  require 
detailed  knowledge  of  subsurface  properties  and  contam¬ 
inant  distributions.  Simplified  screening  models  that  utilize 
upscaled  mass  transfer  coefficients  provide  an  alternative 
method  for  simulating  average  contaminant  behavior  with 
less  site-specific  information  and  computational  cost.  The 
purpose  of  this  work  was  to  evaluate  and  compare  the 
performance  of  such  upscaled  models  in  capturing  average 
concentrations  and  source  longevity  behavior  for  irregular 
source  zones,  and  to  suggest  improvements  that  allow  for 
a  more  general  prediction  capability. 


[36]  Results  from  this  work  demonstrate  that  many  of  the 
existing  simplified  mass  transfer  models  [Dekker,  1996; 
Parker  and  Park ,  2004;  Park  and  Parker,  2005;  Zhu  and 
Sykes,  2004;  Falta  etal.,  2005a]  are  closely  related,  with  the 
main  differences  among  these  models  arising  from  param¬ 
eterization  of  the  upscaled  mass  transfer  correlation.  Com¬ 
parisons  of  predictions  obtained  using  the  upscaled  models 
and  numerical  simulators  show  that  upscaled  models  may 
be  applied  confidently  only  for  the  specific  conditions  under 
which  the  correlation  parameters  were  developed,  and  can 
either  overpredict  or  underpredict  flux-weighted  concentra¬ 
tions  and  source  longevities  by  more  than  one  order  of 
magnitude  when  applied  to  other  scenarios.  This  behavior 
was  primarily  attributed  to  the  effects  of  flow  bypassing  and 
site-specific  saturation  distributions  on  DNAPL  dissolution 
and  mass  recovery. 

[37]  To  address  these  limitations,  existing  upscaled  mod¬ 
els  were  modified  to  replace  the  fitted  upscaled  correlation 
parameters  with  measurable  site-specific  parameters.  Incor¬ 
poration  of  the  initial  average  concentration  and  ganglia- 
to-pool  mass  ratio  into  the  upscaled  mass  transfer  correlation 
was  shown  to  provide  predictions  that  were  similar  (within 
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20%)  to  those  of  the  best  fit  model,  while  maintaining  the 
ability  to  reasonably  predict  (RMSE  <10  mg  L  1 )  flux- 
weighted  concentrations  eluting  from  a  broad  range  of 
source  zone  scenarios.  Limitations  of  the  developed  model 
include  the  overprediction  of  flux-weighted  concentrations 
at  high  levels  of  mass  removal  for  source  zone  scenarios  that 
are  initially  dominated  by  DNAPL  pools,  and  the  need  to 
quantify  the  initial  DNAPL  source  zone  GTP  mass  ratio. 
Overprediction  of  flux-weighted  concentrations  at  high 
levels  of  mass  removal  may  lead  to  underprediction  of 
source  longevity.  Hence,  despite  the  model’s  ability  to 
accurately  predict  flux-weighted  concentrations  at  most 
levels  of  mass  removal  (<85%  mass  removal),  application 
of  the  model  to  predict  posttreatment  source  longevities, 
particularly  at  sites  with  significant  DNAPL  pooling,  should 
be  done  with  caution. 

[3s]  Future  improvement  of  the  predictive  capability  of 
upscaled  mass  transfer  models  will  likely  require  incorpo¬ 
ration  of  additional  site-specific  information  into  the 
upscaled  mass  transfer  correlation  or  modification  of  the 
upscaled  mass  transfer  model  to  account  for  multiple 
domains  possessing  different  mass  transfer  characteristics. 
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